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Abstract 


Euclidean  geometry  is  not  practical  for  describing  complex  natural  shapes.  A 
real-time  solution  using  a  fractal  geometry  approach  is  presented.  First,  existing 
digital  terrain  elevation  data  (DTED)  is  analyzed,  and  then  it  is  interpolated  to 
generate  a  terrain  surface  of  arbitrarily  high  spatial  resolution.  Although  fractal 
geometry  provides  only  an  estimate,  statistical  properties  are  preserved. 
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Introduction 


Fractal  geometry  is  a  mathematical  language  providing  for  concise  descriptions  of  natural 
shapes.  In  particular,  Mark  [1]  has  shown  that  many  geomorphic  surfaces  can  be  modeled 
#  by  a  stochastic  fractal  based  on  fractional  Brownian  motion  (fBm);  fBm  is  a  generalization 

of  classical  Brownian  motion,  where  dilferences  between  successive  positions  are  normally 
distributed.  However,  it  is  necessary  to  restrict  a  single  fBm  function  to  some  range  of  scale 
»  when  simulating  a  topographic  surface. 

This  report  describes  an  application  of  this  theory  to  terrain  modeling.  The  approach 
follows  from  the  work  of  Yokoya  et  al.  [2],  The  next  three  sections  examine  fBm  approxi¬ 
mation,  where  mathematical  purity  is  sacrificed  for  realism  and  execution  speed.  Section  5 
provides  a  specific  example. 

Algorithms  for  the  previous  and  first  level  digital  terrain  elevation  data  (DTED-0)  binary 
data  extraction  [3]  were  written  in  C  and  C++,  respectively.  This  native  code  was  then 
integrated  with  a  graphical  user  interface  (GUI)  and  a  three-dimensional  (3D)  scene  graph 
interactive  display  using  Java  technology.  In  particular,  the  GUI  was  done  using  Java  with 
Swing  components,  and  the  3D  computer  graphics  were  written  in  Java3D.  The  Java  Native 
Interface  (JNI)  allows  for  native  code  to  invoke  the  Java  Virtual  Machine  (JVM)  [4].  Java 
was  chosen  because  of  its  ”  write  once,  run  anywhere”  feature  with  an  acceptable  response 
time. 

2.  Fractional  Brownian  Function 

A  stochastic  function  /(x)  constructs  fractional  Brownian  motion  if,  for  all  x  and  Ax, 

Pr(^y>  <  t)  =  F(t),  (1) 

where  H  is  a  constant  that  lies  in  the  range  0  <  H  <  1  and  F(t)  is  a  cumulative  distribution 
function  of  a  random  variable  t.  Bold-faced  type  indicates  a  vector  quantity;  e.g.,  x  is  the 
position  vector  to  point  x  and  ||  Ax||  is  the  length  of  vector  Ax.  Since  we  are  only  concerned 
with  3D  surfaces,  the  fractal  dimension  of  the  distribution  is  determined  by 

D  =  3  -  H.  (2) 

For  our  particular  application,  equation  (1)  is  the  probability  that  the  difference  in  elevation 
over  some  scaled  distance  is  less  than  t. 

» 

By  assuming  that  F(t)  is  a  zero-mean  Gaussian  distribution  of  variance  a 2,  Yokoya  et 
al.  [2]  rewrites  equation  (1)  as 

l°g  E[\f(x  +  Ax)  -  /(x)|j  =  H  log  ||  Ax||  +  log  C,  (3) 

where  F[|/(x  +  Ax)  —  /(x)|]  is  the  expected  value,  or  simply  the  mean,  of  the  difference  of 
function  values  over  some  distance  ||Ax||,  and  C  =  2a/y/2Tt  (see  [2],  p.  286).  Note  that  for 
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constant  H  and  C,  a  log-log  plot  of  E[\f(x  +  Ax)  —  /(x)|]  as  a  function  of  ||Ax||  has  the 
slope-intercept  form  of  a  straight  line.  This  is  called  a  fractal  plot. 

3.  Fractal  Analysis 

We  cannot  expect  a  geomorphic  surface  to  have  the  characteristics  of  a  true  stochastic 
fractal.  An  exact  fBm  process  is  self-affine,  a  property  that  is  inconsistent  with  most  natural 
surfaces.  In  contrast  to  a  self-similiar  function  where  scale  factors  of  independent  variables 
are  the  same,  self-affine  functions  have  different  scale  factors.  A  natural  fractal,  such  as 
terrain,  should  be  approximated  by  one  or  more  fBm  functions,  each  having  a  range  of  scale. 

Yokoya  describes  a  technique  for  determining  scale  limits  [||Ax||m{n,  ||  Ax||maJ  between 
which  a  natural  fractal  can  be  well-described,  by  a  single  fBm  function.  The  lower  limit 
||  Ax||OT,-n  simply  corresponds  to  the  spatial  resolution  of  the  original  data,  which  is  assumed 
to  be  uniform  and  rectangular.  The  upper  limit  ||Ax||maa.  is  determined  by  a  linearity  test 
of  the  fractal  plot.  For  a  two-dimensional  (2D)  fractal  plot,  the  measure  of  linearity  is  given 
by  the  following  expression: 

J  _  4^11 2+(W0 -w; F  /  ,\ 

—  W0+W2  ’  V  > 

where  /z2 o  is  the  variance  of  log  ||  Ax||,  fi02  is  the  variance  of  log  E[\f(x  +  Ax)  —  /(x)|],  and 
fin  is  the  covariance  of  the  set  of  n  points.  The  upper  limit  of  scale  is  computed  as 

ilAxllma*  =  II  Axllmin +«*-!>  (5) 

where  n*  is  the  minimum  of  n  that  provides  the  local  peak  of  linearity  (see  [2],  pp.  288-291). 

Least-squares  regression  analysis  is  then  used  to  estimate  H,  which  is  the  slope  of  the 
limited  fractal  plot.  The  fractal  dimension  follows  from  equation  (2).  Finally,  the  standard 
deviation  a  =  C >/2x/2  can  be  calculated  since  we  know  log  C,  the  intercept  of  the  least- 
square  line.  The  parameter  H  and  the  standard  deviation  a  of  F ( t )  for  some  range  of  scale 
[|lAxIL»JAxILJ  are  referred  to  as  the  fractal-based  features. 

4.  Fractal  Interpolation 

Assuming  that  statistical  properties  of  a  natural  fractal  are  preserved  at  smaller  scales, 
Yokoya  et  al.[2]  extends  a  range  of  scale  to  [0,  ||Ax||moa,].  His  stochastic  interpolation  al¬ 
gorithm,  called  recursive  midpoint  displacement,  effectively  adds  a  weight  to  the  average  of 
four  neighboring  values  in  computing  the  midpoint.  A  weight  is  the  product  of  extracted 
fractal  features  and  a  random  element. 

Fractal  features  H  and  a  are  determined  by  fractal  analysis  of  the  original  data  (see 
section  3).  It  is  assumed  that  an  N  x  M  rectangular  grid  of  input  data  is  aligned  and  uniform. 
The  values  where  both  i  and  j  are  odd  natural  numbers  satisfying  1  <  i  <  2 N  and 
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Figure  1:  Polygonal  surface. 


1  <  j  <  2 M,  represent  the  known  surface.  First,  compute  values  where  both  i  and  j  are 
even: 

/(bi)  =  4 [/(£  —  l)(i  +  1)  +  f(i  -  1)0’  -  1)  +  f(i  +  1)0  +  1)  +  f(i  +  1)0  _  1)]+  (6) 

[\/l  —  22^/_2  *  ||  AxU^  *  a  *  Gauss ()]. 

When  only  i  or  j  is  even,  the  values  are  calculated  as 

/(m)  =  }[/(*■  - 1)0) + /(oo  + 1)  +  mu  - 1) + /(< + 1)0)]+  (7) 

[2-tf/Vl  -  22^-2  *  ||  Ax||ff  *  a  *  Gauss()]. 

The  function  Gauss()  returns  a  Gaussian  random  number  with  zero  mean  and  unit  variance, 
i.e.,  1V(0, 1).  The  refined  grid  at  this  level  of  recursion  is  (2iV  —  1)  x  (2 M  —  1).  Iteration  of 
this  procedure  results  in  data  of  the  desired  resolution. 

5.  An  Application 

A  wireframe  display  for  a  10  x  5  km  region  of  eastern  USSR  is  shown  in  Figure  1; 
the  precise  location  is  54°0/0”  N  latitude,  36°0/0”  E  longitude.  The  100-m  resolution  grid  of 
DTED  is  aligned  and  uniform.  The  terrain  surface  is  approximated  by  10,000  planar  triangles 
which  were  determined  by  tessellation  of  the  DTED.  *  Intermediate  heights  are  bilinearly 
interpolated  from  neighboring  elevation  posts.  Affine  transformations  of  the  surface  provide 
the  desired  view. 

‘Optimal  triangulation  was  used  to  avoid  generating  triangles  with  long  edges  and  sharp  angles. 
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Figure  2:  Approximated  ffim  surface. 


The  program  provides  for  either  a  wireframe  display  or  Gouraud  solid  area  scan  conver¬ 
sion.  Figure  2  shows  a  shaded  display  of  the  same  region.  In  this  case,  the  original  terrain 
area  was  sampled  every  200  meters  instead  of  every  100  meters.  The  refined  surface  was 
then  generated  from  one  level  of  recursive  midpoint  displacement. 

Although  the  human  eye/brain  combination  is  capable  of  reducing  large  quantities  of  data 
quickly,  it  is  not  as  capable  in  discerning  detail.  A  simple  way  to  compare  the  actual  terrain 
against  the  fractally  generated  terrain  is  a  difference  statistic.  This  statistic  measures  the 
closeness  of  two  datasets  that  are  in  perfect  register  (i.e.,  uniform  cind  aligned  with  elevation 
values  attached  to  the  same  spatial  coordinates).  For  the  previous  example,  the  original 
DTED  has  fj,0  =  233.088  and  a0  =  14.735,  and  the  computed  fractal  terrain  has  fXf  =  233.060 
and  <j j  —  14.568.  This  results  in  a  difference  statistic  with  parameters  [id  =  0.913  and 
ad  =  1.098,  which  suggests  high  correlation. 

Recall  that  the  parameters  H  and  a  characterize  a  fBm  function  ranging  over  [||  Ax||min, 
||  Ax||mar]  (see  section  3).  This  same  function  is  assumed  to  hold  for  the  range  [0,  ||  Ax||min] 
and  allows  for  processing  data  at  higher  levels  of  resolution  if  desired. 

6.  Significance 

Digital  terrain  databases  are  needed  for  analyzing  many  different  battlefield  activities  in 
a  simulation  environment,  e.g.  planning  routes  of  attack  and  determining  line  of  sight  (LOS). 
Many  simulators  use  low  resolution  data  due  to  financial  and/or  hardware  constraints.  This 
limitation  may  be  addressed  through  the  proposed  terrain  synthesis  of  Yokoya  et  al.[2]. 

Our  software  includes  this  procedure  for  real-time  fractal  interpolation  once  two  param¬ 
eters  are  computed  for  a  specific  area  -  the  fractal  dimension  and  standard  deviation.  This 
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approach  also  allows  for  real-time  alteration  of  terrain  within  this  dynamic  environment,  e.g. 
artillery  bombing  which  may  effect  the  soldier’s  decision  making. 

7 .  Future  Work 

More  rigorous  statistical  testing  should  be  considered  when  comparing  actual  data  against 
simulated  terrain  data.  One  possibility  is  a  Kolmogorov- Smirnov  (KS)  test;  KS  is  a  test  of 
distributions  and  thus  more  flexible  than  a  simple  difference  statistic.  Also  being  considered 

is  a  Student’s  t-test,  where  individual  parameters  such  as  the  mean  and  variance  are  more 
closely  examined. 

Computer  ProSram  was  built  on  a  Silicon  Graphics  02  visual  workstation  running 
^  version  °f  the  UNIX  °Perating  system.  Supporting  system  software  includes  a 
compiler/translator  (v2.8.1),  which  is  necessary  for  compiling  native  code, 
and  Java2  Execution  Environment  (Software  Development  Kit  vl.2.1)  with  a  Java3D  API 
(v  1.1.1).  We  are  currently  porting  this  software  to  a  microcomputer  having  a  Win32  GNU 
0/C++  compiler/translator  (available  free  of  charge  from  Cygnus,  Inc.  at  www.cygwin.com) 
and  Java2  1.3,  Java3D  1.2  APIs.  An  effort  is  now  being  made  to  generate  relocatable  dynamic 
link  libraries  (DLLs)  analogous  to  shared  objects  (.so)  in  a  UNIX  environment. 
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